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Abstract. The accuracy of finite-difference time-domain (FDTD) modelling of 
left-handed metamaterials (LHMs) is dramatically improved by using an averaging 
technique along the boundaries of LHM slabs. The material frequency dispersion 
of LHMs is taken into account using auxiliary differential equation (ADE) based 
dispersive FDTD methods. The dispersive FDTD method with averaged permittivity 
along the material boundaries is implemented for a two-dimensional (2-D) transverse 
electric (TE) case. A mismatch between analytical and numerical material parameters 
(e.g. permittivity and permeability) introduced by the time discretisation in FDTD 
is demonstrated. The expression of numerical permittivity is formulated and it is 
suggested to use corrected permittivity in FDTD simulations in order to model LHM 
slabs with their desired parameters. The influence of switching time of source on the 
oscillation of field intensity is analysed. It is shown that there exists an optimum value 
which leads to fast convergence in simulations. 
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1. Introduction 

Recently a great attention has been paid on the research of a new type of artificial 
materials: medium with simultaneously negative permittivity and permeability which 
is introduced by Veselago in his early paper in 1968 [1] and named as left-handed 
metamaterial (LHM). The electric field, magnetic field and wave vector of an 
electromagnetic plane wave in such materials form a left-handed system of vectors. 
The LHMs introduce peculiar yet interesting properties such as negative refraction, 
reversed Doppler effect and reversed Cerenkov radiation etc. One of the most important 
applications of LHMs suggested by Sir John Pendry, is the "perfect lens" [2], e.g. 
the subwavclcngth imaging which allows the information below the diffraction limit 
of conventional imaging systems to be transported. The LHM lenses provide unique 
properties of negative refraction and amplification of evanescent waves, which accounts 
for the reconstruction of source information at the image plane. 

The finite-difference time-domain (FDTD) method [3] is a versatile and robust 
technique. Over the years it has been widely used for the modelling of electromagnetic 
wave interaction with various frequency dispersive and non-dispersive materials. For 
modelling of LHMs with negative material properties, the frequency dispersion has 
to be taken into account therefore the dispersive FDTD method needs to be used. 
The existing frequency dispersive FDTD methods can be categorised into three types: 
the recursive convolution (RC) method [4], the auxiliary differential equation (ADE) 
method [5] and the Z-transform method [6]. The RC scheme relates electric fiux density 
to electric field intensity through a convolution integral, which can be discretised as a 
running sum. The dispersive FDTD method applying the RC scheme has been used 
for modeUing of different types of dispersive materials in [7-14]. The ADE method 
introduces additional differential equations in order to describe frequency dependent 
material properties [15-20]. Another dispersive FDTD method is based on the Z- 
transforms [21,22]: the time-domain convolution integral is reduced to a multiplication 
using the Z-transform, and a recursive relation between electric flux density and electric 
fleld is derived. 

There have been a number of attempts to model LHMs using the FDTD method 
[23-28]. It may seem that the conventional dispersive FDTD has been verified in the 
literature: the negative refraction effect which is inherent to the boundary between the 
free space and LHM is observed and the planar superlens behaviour has been successfully 
demonstrated [23-25] . Actually, this means that the LHM is correctly modelled only for 
the case of propagating waves. When evanescent waves are considered the conventional 
implementation of dispersive FDTD method may lead to inaccurate results. Usually, the 
evanescent waves decay exponentially over distances and thus they are concentrated in 
the close vicinity of sources, that is why conventional FDTD modelling of non-dispersive 
materials does not suffer from the aforementioned numerical inaccuracy. In the case of 
LHM, the evanescent waves play a key role and have to be modelled accurately because 
of the perfect lens effect [2]. This explains why early FDTD simulations have not 
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demonstrated the subwavelength imaging property of LHM lenses [23,24]. A slab of 
LHM effectively amplifies evanescent waves which normally decay in usual materials 
and allows transmission of subwavelength details of sources to significant distances. 

Other numerical studies using the FDTD method include the effect of losses and 
thicknesses on the transmission characteristics of LHM slabs [27], and the influence 
of numerical material parameters on their imaging properties [28] etc. Besides the 
FDTD method, the pseudo-spectral time-domain (PSTD) method has been used for the 
modelling of backward- wave metamaterials [29]. It is claimed in [29] that the FDTD 
method cannot be used to accurately model LHMs due to the numerical artefact of 
the staggered grid in FDTD domain. However, we shall show later by comparing the 
transmission coefficient calculated from FDTD simulation and exact analytical solutions 
that with proper field averaging techniques [28,30], the FDTD method indeed can be 
used to accurately characterise the behavior of both propagating and evanescent waves 
in LHM slabs. Furthermore, it has been reported in [31,32] that with special treatment 
(i.e. averaging techniques) along material boundaries, accurate modelling of curved 
surfaces of conventional dielectrics as well as surface plasmon polaritons between metal- 
dielectric interfaces can be achieved without using extremely fine FDTD meshes. 

Ideally lossless LHM slabs with infinite transverse length provide unhmited 
subwavelength resolution. However in realistic situations, the subwavelength resolution 
of the LHM lenses is limited by losses [33], the thickness of the slab and the mismatch 
of the slab with its surrounding medium [34]. It is important to understand these 
theoretical limitations because they can help verify numerical simulations. In this 
paper, we have performed the modelling of infinite LHM slabs and their transmission 
characteristics. The infinite LHM slab is modelled using the periodic boundary condition 
and a material parameter averaging technique is used along the boundaries of LHM slabs. 
In contrast to FDTD modelling of conventional dielectric slabs where the averaging is 
only a second-order correction to improve the accuracy of simulations, the averaging of 
permittivity is an essential modification for modelling of LHM slabs. The averaging 
of material parameters implemented in our FDTD simulations is equivalent to the 
averaging of current density originally introduced in [28] and is analysed in detail in 
this paper. It is demonstrated that other numerical aspects such as numerical material 
parameters and the switching time of source also have considerable infiuences on FDTD 
simulations. 

2. Dispersive FDTD Modelling of LHMs with Spatial Averaging at the 
Boundeiries 

We consider here lossy isotropic LHM slabs modelled using the effective medium method. 
The Drude model is used for both the permittivity ^(a;) and permeability /^(a;) with 
identical dispersion forms: 
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/i(c^) = /io f 1 - , ) , (2) 

where cOpe and cUpm are electric and magnetic plasma frequencies and 7e and 7m are 
electric and magnetic collision frequencies, respectively. 

Although there are various dispersive FDTD methods available for the modelling 
of LHMs, due to its simplicity and efficiency, we have implemented the ADE method 
in this paper. There are also different schemes involving different auxiliary differential 
equations in addition to conventional FDTD updating equations. In this paper, two 
schemes, namely the (E, J, H, M) scheme [3] and the (E, D, H, B) scheme [5], are 
used and introduced respectively. 

2.1. The (E, D, H, B) Scheme 

The (E, D, H, B) scheme is based on Faraday's and Ampere's Laws: 

5B 

curl(E) = - ^, (3) 

curl(H) = — , (4) 

as well as the constitutive relations D = eE and B = /xH where e and are expressed 
by ([I]) and ([2]), respectively. Equations ([3]) and (jl]) can be discretised following a normal 
procedure [3] which leads to conventional FDTD updating equations: 

gn+i = B" - At ■ ^(E"+5), (5) 
j^n+l = 13" + ^^. ^(H"+^ ) . (6) 



where curl is discrete curl operator. At is FDTD time step and n is the number of time 
steps. 

In addition, auxiliary differential equations have to be taken into account and they 
can be discretised through the following steps. The constitutive relation between D and 
E reads 

(u^ - JU-fe) D = £o - J^le ' ^'e) E. (7) 

Using inverse Fourier transform and the following rules: 

9 2 9^ 
^ 5^' " -9^' 
Equation (171) can be rewritten in the time domain as 

The FDTD simulation domain is represented by an equally spaced three- 
dimensional (3-D) grid with periods Ax, Ay and Az along x-, y- and 2;-directions, 
respectively. For discretisation of ([9]), we use central finite difference operators in time 
{St and Sf) and central average operator with respect to time {^t and fif): 
92 52 d St 222 



dt^ ^ (At)2' dt ^ At^'' "^P' ^ "^f^^*' 



where the operators 5i, 5^, fit and /i^ are defined as in [35]: 
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Here F represents field components and m^;, m^, are indices corresponding to a certain 
discretisation point in FDTD domain. The discretised Eq. ^ reads 
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Note that in (ITTl) . the discretisation of term of is performed using the central 
average operator /i^ in order to guarantee improved stability; the central average 
operator /if is used for the term containing 7e to preserve second-order feature of the 
equation. Equation ffTTj) can be written as 
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Therefore the updating equation for E in terms of E and D at previous time steps is as 
follows: 
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The updating equation for H is in the same form as f[T^ by replacing E, D, u. 
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Equations ([5]), ([6]), (fT3|) and (fT4|) form an FDTD updating equation set for LHMs using 
the (E, D, H, B) scheme. If both the plasma frequency and collision frequency are 
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equal to zero i.e. ujpe = ujpm = and 7e = 7m = 0, then they reduce to the updating 
equations in the free space. 



2.2. The (E, J, H, M) Scheme 

An alternative ADE FDTD scheme starts with different forms of Faraday's and Ampere's 
Laws for LHMs: 
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Following the same procedure as for the (E, D, H, B) scheme, Eqs. f|T5l) - f|T8l) can be 
discretised as: 

At 



(9E 

curl(H) = sq— + J, 
where the electric and magnetic current density, J and M are defined as 
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Again Eqs. (fT9l) - (!22l) become the free space updating equations if both the plasma 
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2.3. The Spatial Averaging Methods 

In addition to the above introduced ADE schemes, due to the staggered grid in 
FDTD domain, a modification at the interfaces between different materials is often 
used to improve the accuracy of FDTD simulations. It has been shown that 
the field averaging techniques based on the averaging of material parameters (e.g. 
permittivity and permeability) provide a second-order accuracy [36]. The averaged 
permittivity/permeability can be obtained by performing either arithmetic mean, 
harmonic mean or geometrical mean [36] and the arithmetic mean has been proven to 
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Figure 1. The layout of FDTD grid illustrating the arrangement of material 
boundaries along j/-direction. The gray arrows indicate where the averaged permittivity 
is used. The FDTD unit cell is shown on the right side. 



have the best performance amount these three schemes. Previous analysis of averaging 
techniques are performed for conventional dielectrics with positive permittivity and 
permeability. For the materials with negative permittivity/permeability, one of the 
simplest ways to implement averaging is to use arithmetic mean. Furthermore, averaging 
should be applied only for the field components tangential to material interfaces. 
Therefore depending on the configuration of FDTD simulation domain e.g. two- 
dimensional (2-D) TE, 2-D TM or three-dimensional (3-D) cases, the averaging needs to 
be performed in different ways. In this paper we have considered a 2-D (x-y) simulation 
domain with H-polarisation where H is directed only along 2;-direction. Therefore only 
three field components are non-zero: Er^, Ey and Hz- For the interfaces between LHM 
slab and the free space along x-direction, the averaged permittivity for the tangential 
electric field component E^ is given by 
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which is equivalent to replacing the plasma frequency Upe by u'^^ 
Therefore along the boundaries, the updating equation for E^ reads 
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The locations where the updating equation fl2^ is used is illustrated in Fig. [Hby gray 
arrows. 

The averaging of permittivity can be implemented for the (E, D, H, B) scheme. 
While for the (E, J, H, M) scheme, it is proposed in [28] to use the averaging of 
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tangential current density along the boundaries of LHM slab. The averaged current 
density can be calculated as (the free space current density Jq = 0): 

< J. >= = y, (25) 

then the updating equation for along the boundaries of LHM slab becomes (expanded 
from Eq. (|20l)) 

^x\mx,my + A {^Hz\m^,my ~ I .m^ - 1 ) (26) 
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Theoretically the above two averaging methods have same effects due to the linear 
relations 

B = eE = eoE + —J, B = /iH = /ioH + — M. (27) 
jtc; juj 

Therefore the averaging of current density is identical to the averaging of permeability. 

In this paper, we have used the (E, D, H, B) scheme in all our simulations because of 

its simplicity in implementation. In order to demonstrate the advantage of averaging 

technique, we have also compared the results from simulations with and without 

averaged permittivity along material boundaries. For the case of without averaging, 

the tangential electric fields indicated by gray arrows in Fig. [1] are updated using their 

updating equations in the free space. 

The above averaging of permittivity only applies to the field components tangential 

to material interfaces and for the case of TE polarisation considered in our simulations. 

If it is required to apply the averaging schemes to materials with planar boundaries for 

TM and three-dimensional (3-D) cases or even for structures with curved surfaces, one 

can follow the procedures introduced in [31,32]. 



3. Numerical Implementation 

For simplicity, in our simulations we assume that the plasma frequency is Upe = ujpm = 
Up = \/2caj where to is the operating frequency, therefore matched LHM slabs are 
modelled in our simulations. A small amount of losses is used i.e. 7e = 7m = 7 = 
0.0005a; which gives relative permittivity and permeability = /^r = — 1 — O.OOlj to 
ensure the convergence of simulations. It is worth mentioning that there is a small 
amount of mismatch between numerical (in FDTD domain) and analytical permittivity 
([1]) which is caused by FDTD time discretisation [28] . However, such a mismatch causes 
the amplification of transmission coefficient only for lossless LHM slabs or when the 
losses are very small. For the amount of losses used in our simulations, the effect of 
mismatch is damped and no amplification is found in transmission coefficient. The effect 
of FDTD cell size on this mismatch is analysed in later sections. 
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Figure 2. Schematic diagram of two-dimensional (2-D) FDTD simulation domain for 
calculation of numerical transmission coefficient 



As shown in Fig. [21 an infinite LHM slab is modelled by applying the Bloch's 
periodic boundary conditions (PBCs). For any periodic structures, the field at any time 
satisfies the Bloch theory, i.e. 

E(x + L) = E(x)e^'^^^, H(a; + L) = H(a;)e^'^^^, (28) 

where x is any location in the computation domain, is the wave number in x-direction 
and L is the lattice period along the direction of periodicity. When updating the fields 
at the boundary of computation domain using the FDTD method, the required fields 
outside the computation domain can be calculated using known field values inside the 
domain through ( l28l) . Since infinite structures can be truncated with any period, for 
saving computation time, we have used only four FDTD cells in x-direction (L = 4Ax). 
Along y-direction, the Berenger's original perfectly matched layer (PML) [37] is used 
for absorbing propagating waves {k^ < /cq), and the modified PML [38] is used when 
calculating the transmission coefficient for evanescent waves [k^ > k^). A soft plane- 
wave sinusoidal source (which allows scattered waves to pass through) with phase delay 
corresponding to different wave number is used for excitations, 

= {t,js) + s{t)e-^'^^''\ (29) 

where js is the location of source along y-direction, s{t) is a time domain sinusoidal 
wave function, i G [1, /] is the index of cell location and / is the total number of cells in 
x-direction (J = 4 in our case). By changing the values of wave number kx-, either pure 
propagating waves [k^ < k^) or pure evanescent waves {k^ > fco) can be excited. 

The spatial resolution in FDTD simulations is Ax = Ay = A/100 where A is the free 
space wavelength at the operating frequency. According to the stability criterion [3] , the 
discretised time step is At = Ax/ ^/2c where c is the speed of light in the free space. As 
illustrated in Fig. [2], the source plane is located at a distance of d/2 to the front interface 
of the LHM slab where d is the thickness of the slab. Therefore the first image plane 
is at the centre of the LHM slab and the second image plane is at the same distance 
of d/2 beyond the slab. The spatial transmission coefficient is calculated as a ratio of 
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Figure 3. Comparison of transmission coefficient of infinite planar LHM slabs 
calculated from exact analytical solutions and dispersive FDTD method with and 
without averaging of permittivity along the boundaries of LHM slabs. 



the field intensity at the second image plane to the source plane for different transverse 
wave numbers kx after the steady-state is reached in simulations. 

Figure [3] shows the transmission coefficient for an infinite planar LHM slab with 
thickness d = 0.2A calculated using the FDTD method with and without averaging of 
permittivity along the boundaries, and its comparison with exact analytical solutions. 
It can be seen that using the arithmetic mean of permittivity, the numerical results show 
excellent agreement with analytical solution using spatial resolution Ax = A/100. Good 
correspondence can be also obtained when the FDTD cell size is increased to As = A/80. 
On the other hand, without averaging, the material boundary is not correctly modelled 
which introduces an amplification (resonance) at a location of approximately = 2AkQ 
in transmission coefficient for the case of Ax = A/100. Reducing FDTD cell size to 
Ax = A/200 and Ax = A/400, the behaviour of the resonance remains similar but the 
location shifts to k^ = 2.8/co and k^ = 3.2fco, respectively. Therefore we predict that 
only if a very small FDTD cell size is used in simulations that the results can converge 
to the right solution. Such a comparison demonstrates the significance of the averaging 
technique. Conventionally the arithmetic averaging is only a second-order correction 
for modelling of conventional dielectric slabs, however it is shown in Fig. [3] that for 
modelling of LHM slabs, the averaging becomes an essential modification. 

The results shown in Fig. [3] may explain some incorrect results obtained previously. 
For instance, the amplification of transmission coefficient in [26,27] is caused by incorrect 
modelling of material boundaries, but such amplification is pure numerical and does not 
exist in actual LHM slabs [33,34]. It is claimed in [39] that the imaging property of 
finite-sized LHM slabs is significantly affected by their transverse dimensions, which we 
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suppose that the conclusion is drawn from incorrect numerical simulations. We have 
performed accurate simulations using averaging of material properties and confirmed 
that the resolution of a near-field lens using LHMs is free from its transverse aperture 
size [30]. 

In our simulations for the calculation of transmission coefficient, we have used 
PBCs in x-direction to model infinite structures and averaged permittivity along the 
boundaries in y-direction. If one requires to model finite-sized structures (in both 
X- and y-directions) , the averaged permittivity/permeability needs to be used for the 
corresponding tangential component along the boundaries in both directions. 

Besides the averaging technique used along the boundaries of LHM slabs, there are 
other numerical aspects in FDTD simulations in order to model the behaviour of LHM 
slab correctly and accurately. These aspects are introduced respectively in following 
sections. 

4. Effects of Numerical Material Parameters 

Usually for modelling of conventional dielectrics, the results are assumed to be accurate 
enough i.e. the effect of numerical material parameters can be ignored if an FDTD cell 
size of smaller than Aa; = A/20 is used. However since the discretisation introduces a 
mismatch between numerical and analytical permittivity/permeability, when modelling 
LHMs, especially when the evanescent waves are involved, the FDTD spatial resolution 
has a significant impact on the accuracy of simulation results. The effect of numerical 
permittivity/permeability is originally reported in [28] for lossless LHMs using the (E, 
J, H, M) scheme. Following the same procedure one can also obtain the numerical 
permittivity/permeability for the case of lossy LHMs. In this paper, the numerical 
permittivity/permeability (for lossy LHMs) is derived for the (E, D, H, B) scheme. 
In the case of plane waves, when 

E" = Ee^""'^*, D" = De^"'^^*, (30) 

Eq. (fT3ll reduces to D" = eE", where e is the numerical permittivity of the following 
form: 

a;2(At)2cos2^ 1 , ^ 

[ 2sin^ (2sin^^ -j7Atcos^if^)J 

If the collision frequency 7 = 0, then (|3T|) reduces to the numerical permittivity for 
lossless LHMs given in [28]. 

Previously we have used an FDTD cell size of Ax = A/100 in simulations. 
Substitute the corresponding time step At = Aa:/\/2c and the operating frequency, we 
can obtain the numerical relative permittivity from fl3T|) as e^- = —0.9993 — O.OOlOj. 
Although there is a mismatch between the real part of relative permittivity and 
— 1, the loss in LHMs damps such a mismatch and the simulation results show very 
good accuracy. However if we increase FDTD cell size, the numerical permittivity 
introduces severer mismatch which causes the discrepancy between FDTD simulation 
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Figure 4. Transmission coefficient of infinite planar LHM slabs using the proposed 
FDTD method with averaged permittivity and without the correction of material 
parameters. The amplification of transmission coefficient is caused by the mismatch 
introduced by the time discretisation in FDTD (e,. = -0.9959 - O.OOlOj) with 
Aa; = A/40. The same permittivity is used to obtain the analytical solution for 
comparison. 



result and exact solutions. For example, for the case of Aa; = A/40, the mismatch 
brings an amplification in transmission coefficient as shown in Fig. HI Again using 
( !3T|) we can estimate this mismatch and the numerical relative permittivity reads 
Sr = —0.9959 — O.OOlOj. Using such permittivity in analytical formulations, we can 
obtain the corresponding transmission coefficient which is also plotted in Fig. H] for 
comparison. A good correspondence is shown and at high-wave-vector region, the 
discrepancy is cause by insufficient sampling points as for the case of using large cell 
size (e.g. Ax > A/10) for conventional FDTD. 

Another advantage of estimating the numerical permittivity is the correction of 
the mismatch for FDTD simulations. After simple derivations, we can obtain corrected 
plasma frequency and collision frequency as 



„ 2 sin 

— 



-2(e; - 1) sin ^ - e';-fAt cos 



2 



P (At)2 C0S2 ^ 



^ ~ (4-l)Atcos^' ^^^^ 

where and e'^ are the real and imaginary parts of the design relative permittivity Sr, 
respectively. For the case of Er = —I — O.OOlj, substitute = —1 and e" = —0.001 
into (|32|) we get Up = 1.41570; and 7 = 5.0051 x lO"^^;. Using the corrected material 
parameters, the FDTD simulation result and its comparison with analytical solutions 
are shown in Fig. [5l It can be seen that the mismatch has been canceled in FDTD 



13 



(U 

o 
o 

a 
o 



X FDTD 
Analytical 




Figure 5. Transmission coefficient of infinite planar LHM slabs using the proposed 
FDTD method with averaged permittivity and with the correction of material 
parameters. The numerical permittivity in FDTD (Ax = A/40) is = — 1 — O.OOlj. 
The same permittivity is used to obtain the analytical solution for comparison. 



simulations hence there is no amphfication in transmission coefficient. Again the 
discrepancy with exact solutions in high-wave-vector region is caused by insufficient 
sampling points for such an FDTD spatial resolution of Ax = A/40. Therefore we 
suggest to use FDTD cell size smaller than Ax = A/80 for modelling of LHMs especially 
when evanescent waves are involved. 

5. Effects of Switching Time 

Conventionally for single frequency simulations, the source should be smoothly switched 
to its maximum value in order to avoid exciting other frequency components [23]. 
For modelling of LHMs, the switching time has even more significant effect on their 
behaviour. It is well known that the switching time considerably influences the 
oscillation of images and often thirty period is used as the switching time [25, 27]. 
However, perhaps this is the reason that no stable images could be obtained in [23, 27] 
since recently, it is reported in [40] that using a switching time of at least one hundred 
periods one can obtain stabilised image for lossless LHMs. 

In our FDTD simulations, we also notice that switching time influences the 
oscillation of the field intensity at the image plane and hence the convergence time 
in simulations. We have performed FDTD simulations with different switching time 
equal to SOTq, ISOTq and 250To where Tq is the period of the sinusoidal signal. The 
FDTD cell size is Ax = A/100 and corrected material parameters from (l32l) are used. 
In order to ensure faster convergence, we have chosen larger amount of losses and used 
ir = —1 — O.Olj in simulations. It should be noted that because high- wave- vector 
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Figure 6. The influence of different switching time on the convergence time in FDTD 
simulation of infinite LHM slabs for a fixed wave number = Sfco- The Tq is the 
period of the sinusoidal wave function at the operating frequency. The field intensity 
is taken at the second image plane of the LHM slab. 

components travel very slowly in LHM slabs and the process of the growth of evanescent 
waves requires a very long time to reach the steady-state, field values should be taken 
only after total convergence is reached in simulations. For our case of ir = —1 — O.Olj, we 
have used a criteria of 0.001% for detecting iteration errors and terminating simulations. 
It is clearly shown in Fig. [6] that for a fixed wave number {k^ = 3A;o); the oscillation of 
field intensity can be significantly suppressed by prolonging the switching time. 

It is understandable that when the oscillation can be neglected, the convergence 
time increases with the switching time. For demonstration of the impact of switching 
time on convergence time, we have performed FDTD simulations with various switching 
time. The collected data is plotted in Fig. [7l It can be seen that there exists an optimum 
switching time when the minimum convergence time can be achieved for the case of 
kx = 3/co- However for different wave vectors and different material parameters, the 
behaviour of oscillation differs considerably and in certain cases the oscillation may last 
until a very long time. For practical simulations such as modelling of subwavelength 
imaging by a line source, since the source contains all wave vectors, therefore it is 
necessary to switch the source slowly enough to ensure and speed up the convergence 
of simulations. 

6. Conclusions 

In conclusion, we have performed simulations of LHMs using the dispersive FDTD 
method. Two ADE methods namely the (E, D, H, B) scheme and the (E, J, H, M) 
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Figure 7. The dependence of the convergence time on the switching time in FDTD 
simulations of infinite LHM slabs for a fixed wave number = 3fco. A criteria of 
0.001% is used to detect iteration errors and terminating simulations. 



scheme which lead to exactly same results and the respective averaging techniques along 
material boundaries are introduced. The comparison with exact analytical solutions 
demonstrates that the averaging of permittivity/permeability along the boundaries of 
LHM slabs is essential for correct and accurate modelling of LHMs. The numerical 
permittivity in FDTD is formulated where a mismatch between numerical and analytical 
permittivity is introduced by FDTD time discretisation. We suggest to correct such 
a mismatch in order to model LHMs with their desired parameters in FDTD. The 
behaviour of oscillation of field intensity for different switching time is also analysed. 
It is shown that there exists an optimum value which leads to fast convergence in 
simulations. 
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